setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/varscan/bayes_p_errobar")
sel=list.files(pattern = ".txt")
library(ggplot2)
result=as.data.frame(matrix(NA,1,2))
names(result)=c("unitID","group")
result=result[-1,]

for(i in 1:length(sel)){
  file=read.table(sel[i],head=T,sep="\t")
  file1=tidyr::separate(file,normal_rt,into=c("normal_bayes_pvalue","normal_beta0","normal_lower","normal_upper"),sep="_")
  file2=tidyr::separate(file1,tumor_rt,into=c("tumor_bayes_pvalue","tumor_beta0","tumor_lower","tumor_upper"),sep="_")
  file2$total_reads=as.numeric(file2$normal_reads1)+as.numeric(file2$normal_reads2)+as.numeric(file2$tumor_reads1)+as.numeric(file2$tumor_reads2)
  file3=file2[file2$total_reads>=10,]
  file3$normal_bayes_pvalue=as.numeric(file3$normal_bayes_pvalue)
  file3$tumor_bayes_pvalue=as.numeric(file3$tumor_bayes_pvalue)
  file3$unitID=paste(file3$chrom,file3$position,sep=":")
  bias=file3[file3$normal_bayes_pvalue<0.05|file3$tumor_bayes_pvalue<0.05,]
  bias_normal=data.frame(VAF=bias$normal_var_freq,reads1=bias$normal_reads1,reads2=bias$normal_reads2,pvalue=bias$normal_bayes_pvalue,beta=bias$normal_beta0,lower=bias$normal_lower,upper=bias$normal_upper,unitID=bias$unitID)
  bias_tumor=data.frame(VAF=bias$tumor_var_freq,reads1=bias$tumor_reads1,reads2=bias$tumor_reads2,pvalue=bias$tumor_bayes_pvalue,beta=bias$tumor_beta0,lower=bias$tumor_lower,upper=bias$tumor_upper,unitID=bias$unitID)
  bias_normal=bias_normal[bias_normal$pvalue<0.05,]
  bias_tumor=bias_tumor[bias_tumor$pvalue<0.05,]
  biasdata=rbind(bias_normal,bias_tumor)
  biasdata$VAF=as.numeric(gsub("%","",biasdata$VAF))/100
  biasdata$group=ifelse(biasdata$reads1>biasdata$reads2,"down","up")
  rt=data.frame(unitID=biasdata$unitID,group=biasdata$group)
  result=rbind(result,rt)
  
}

id=unique(result$unitID)
frt=as.data.frame(matrix(NA,1,3))
names(frt)=c("unitID","down","up")
frt=frt[-1,]

for(j in id){
tp=result[result$unitID==j,]
tmp=data.frame(unitID=j,down=table(tp[tp$group=="down",])[1],up=table(tp[tp$group=="up",])[1])
frt=rbind(frt,tmp)
}
frt[is.na(frt)]=0
frt$group=ifelse(frt$down>frt$up,"down","up")
write.csv(frt,"E:/5hmc_file/H3k的分析/bias_up_down_statis.csv",quote=F,row.names = F)







######################################################################
setwd("E:/5hmc_file/H3k的分析")
deepsea=read.table("tmp/d85a2d61-a41f-4e25-8199-b3a319ea96c3_727vcf_FEATURE_zscore.tsv",head=T,sep="\t")
coln=names(deepsea)
AShm=read.csv("bias_up_down_statis.csv",head=T)
AShm=AShm[,-c(2:3)]

braincell=c("Brain_Angular_Gyrus","Brain_Anterior_Caudate","Brain_Germinal_Matrix","Brain_Hippocampus_Middle","Brain_Inferior_Temporal_Lobe","Brain_Mid_Frontal_Lobe","Brain_Substantia_Nigra","Fetal_Brain_Male","Fetal_Brain_Female")

for(i in 1:length(braincell)){
  braincn=coln[grep(coln,pattern = braincell[i])]
  braincn=braincn[grep(braincn,pattern = "H3K")]
  deepseq=deepsea[,c("chrom","pos",braincn)]
  result=data.frame(matrix(NA,727,length(braincn)))
  H3k=unlist(strsplit(braincn,"\\."))[seq(2,3*length(braincn),3)]
  colnames(result)=H3k
  
  frt=data.frame(matrix(NA,727,length(braincn)))
  colnames(frt)=H3k
  
  
  deepseq$unitID=paste0(deepseq$chrom,":",deepseq$pos)
  deepseq=merge(deepseq,AShm,by="unitID")
  for(j in 1:length(braincn)){
    result[,j]=ifelse(deepseq[,3+j]<0,"down","up")
    frt[,j]=ifelse(result[,j]==deepseq$group,"same","opposite")
  }
  
  rt_statis=data.frame(matrix(NA,length(braincn),4))
  row.names(rt_statis)=H3k
  colnames(rt_statis)=c("same","opposite","same_ratio","pvalue")
  for(k in 1:length(braincn)){
    rt_statis[k,1]=table(frt[,k]=="same")[2]
    rt_statis[k,2]=table(frt[,k]=="opposite")[2]
    rt_statis[k,3]=rt_statis[k,1]/(rt_statis[k,1]+rt_statis[k,2])
    rt_statis[k,4]=binom.test(rt_statis[k,1],rt_statis[k,1]+rt_statis[k,2],p=0.5)$p.value
  }
  fn=paste0(braincell[i],"_statis.csv")
  write.csv(rt_statis,fn,quote=F,row.names = T)
}
